//========================================================================================
// Sequence3Bit
// ----------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/src/Sequence3Bit.cc#6 $
// ------------------------------------------------------------------------------
// $Log: $
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//


#include "Sequence3Bit.hh"


void 
Sequence3Bit::Bug(int n)
{
  cout<<"You've found a bug n="<<n<<" mBaseSize="<<mBaseSize<<endl;
}

//----------------------------------------
// * Init
//----------------------------------------
//
void
Sequence3Bit::Init(){
  mBaseSize=0;mMerIdx = 0;mBaseIdx=0;

  // Precompute some of the values we need to reduce the 
  // number of operations required on the fly.  Surprisingly, 
  // This doesn't actually seem to save any time.  I would 
  // have thought that 3 array lookups would be much quicker
  // than 1 divide, two multiplies, one subtract, and a shift, 
  // 5 math operations, vs 3 array lookups, but it's a wash. 
  // If array lookups were worse, I'd thiunk it was a cache
  // effect, math is all in cache, array lookups are not. 
  // 
  for(int i = 0;i< PRECOMPUTE_TABLE_SIZE;i++){
    //	 Sequence3Bit::WordNum[i] = i /10;
    //Sequence3Bit::shiftAmt[i] = (i - (Sequence3Bit::WordNum[i] * 10))*3;
    //Sequence3Bit::Mask3Bit[i] = 0x07 << Sequence3Bit::shiftAmt[i];
  }
}

int 
Sequence3Bit::NumBases(){ return(mBaseSize);};


//----------------------------------------
// * OnlyBases
//----------------------------------------
// Returns true if the sequence in the range [start,start+chkSize-1] contains only 
// bases.  Assumes sequence is element of {A,C,T,G,N,X,a,c,t,g,n,x}.  Masks the ASCII
// to determine if each character is a base rather than N,X,n,x.  Bit 3 is diagnostic
// as follows:
//
// a = 61 (hex)    97 (dec)      0110 0001   0110 0(00)1   a -> 0
// c = 63 (hex)    99 (dec)      0110 0011   0110 0(01)1   c -> 1
// g = 67 (hex)    103 (dec)     0110 0111   0110 0(11)1   g -> 3
// t = 74 (hex)    116 (dec)     0111 0100   0111 0(10)0   t -> 2
// A = 41 (hex)    65 (dec)      0100 0001   0100 0(00)1   A -> 0
// C = 43 (hex)    67 (dec)      0100 0011   0100 0(01)1   C -> 1
// G = 47 (hex)    71 (dec)      0100 0111   0100 0(11)1   G -> 3 
// T = 54 (hex)    84 (dec)      0101 0100   0101 0(10)0   T -> 2 
// n = 6e (hex)    110 (dec)     0110 1110   0110 (111)0   n -> 7
// x = 78 (hex)    120 (dec)     0111 1000   0111 (100)0   x -> 4
// N = 4e (hex)    78 (dec)      0100 1110   0100 (111)0   N -> 7
// X = 58 (hex)    88 (dec)      0101 1000   0101 (100)0   X -> 4  
//
bool 
Sequence3Bit::OnlyBases(int start,int chkSize){
  register int end = start + chkSize - 1;
  if (end >= mBaseSize) end = (mBaseSize-1);

  // In 3Bit encoding X and N are 1 at bit position 2 in 3 bit encoding. 
  // NORX3BITMASK = 0x04
  for(register int i = start;i<= end;i++){
    if ((*this)[i] & NORX3BITMASK) return(false);
  }
  return(true);
}

//----------------------------------------
// * ContainsNX
//----------------------------------------
// Returns true if the sequence in the specified range contain any N's or X's. 
// 
bool 
Sequence3Bit::ContainsNX(int start,int chkSize){
  register int end = start + chkSize - 1;
  if (end >= mBaseSize) end = (mBaseSize-1);

  // In 3Bit encoding X and N are 1 at bit position 2
  for(register int i = start;i<= end;i++){
    if ((*this)[i] & NORX3BITMASK) return(true);
  }
  return(false);
}

//----------------------------------------
// * Print
//----------------------------------------
// Prints the sequence to the given stream (which can be cout).  
// 
void 
Sequence3Bit::Print(ostream &fout){

  if (mBaseSize ==0){
    cout<<"No sequence"<<endl;
    return;
  }

  for(int i =0;i< mBaseSize;i++){
    unsigned int baseNum = (*this)[i];
    char c = Bit2Base(baseNum);
    fout<<c;
    if (!(i%50) && (i!= 0)) cout<<endl;
  }
}


//----------------------------------------
// * Set
//----------------------------------------
// Sets base specified by (wordIdx,bitIdx) to the base given in ASCII by c.  
// wordIdx is the int within the Sequence3Bit that contains the base, and 
// bitIdx is which set of 3-bits within that word.  Each int, you will recall, stores
// ten bases at 3-bits/base. (30 of 32 bits used).  When reading in a series
// of bases, it is easy to increment wordIdx and bitIdx, so that is why 
// set takes these as input rather than the global baseIdx. 
//
// a = 61 (hex)    97 (dec)      0110 0001   0110 0(00)1   a -> 0
// c = 63 (hex)    99 (dec)      0110 0011   0110 0(01)1   c -> 1
// g = 67 (hex)    103 (dec)     0110 0111   0110 0(11)1   g -> 3
// t = 74 (hex)    116 (dec)     0111 0100   0111 0(10)0   t -> 2
// A = 41 (hex)    65 (dec)      0100 0001   0100 0(00)1   A -> 0
// C = 43 (hex)    67 (dec)      0100 0011   0100 0(01)1   C -> 1
// G = 47 (hex)    71 (dec)      0100 0111   0100 0(11)1   G -> 3 
// T = 54 (hex)    84 (dec)      0101 0100   0101 0(10)0   T -> 2 
// n = 6e (hex)    110 (dec)     0110 1110   0110 (111)0   n -> 7
// x = 78 (hex)    120 (dec)     0111 1000   0111 (100)0   x -> 4
// N = 4e (hex)    78 (dec)      0100 1110   0100 (111)0   N -> 7
// X = 58 (hex)    88 (dec)      0101 1000   0101 (100)0   X -> 4  
//
void 
Sequence3Bit::Set(int wordIdx,int bitIdx, char c){

  // KJD Horrible kludge.  Not sure why my previous resize math is sometimes off. 
  if (mData.size() == wordIdx) mData.resize(mData.size()+1);

  unsigned int val = c;
  register unsigned int word = mData[wordIdx];
  
  // Clear whatever was already there
  // #define LOW3BITSMASK 0x00000007
  register unsigned int mask = LOW3BITSMASK;
  unsigned int bitShift = bitIdx * 3;
  mask = mask << bitShift;
  mask = ~mask;
  word = word & mask;
  
  // Extract encoding from ASCII
  // ASCII3BITMASK 0x0E
  register unsigned int base = (c & ASCII3BITMASK) >> 1;   
  base = base << bitShift;
  
  word = word | base;
  mData[wordIdx] = word;
}

//---------------------------------------------------------------------------
// * Reverse
//---------------------------------------------------------------------------
// Sets Rev to the reverse of (*this).   
// KJD Must write faster version of this.  This is very QND. 
//
void
Sequence3Bit::Reverse(Sequence3Bit &Rev)
{  
  Sequence temp;
  temp.resize(mBaseSize);
  for(int i = 0;i< mBaseSize;i++){
    temp[i] = Bit2Base((*this)[i]);
  }
  temp.Reverse();
  Rev.clear();
  Rev.ConcatenateSequence(temp);
}

//---------------------------------------------------------------------------
// * ReverseComplement
//---------------------------------------------------------------------------
// Sets Rev to the reverse complement of (*this).   
// KJD Must write faster version of this.  This is very QND. 
//
void
Sequence3Bit::ReverseComplement(Sequence3Bit &Rev)
{
  Sequence temp;
  temp.resize(mBaseSize);
  for(int i = 0;i< mBaseSize;i++){
    temp[i] = Bit2Base((*this)[i]);
  }
  temp.ReverseComplement();
  Rev.ConcatenateSequence(temp);
}

//----------------------------------------
// * ConcatenateCString
//----------------------------------------
// Adds a C-language string of ASCII bases to the Sequence3Bit.   
// 
void 
Sequence3Bit::ConcatenateCString(const char *cstring){
  register int clen = strlen(cstring);
  register int oldMerSize = mData.size();
  register long idx;
  
  // Make space for it.
  int tailAmt = 10 - mMerIdx - 1;  
  int newMers = ((clen-tailAmt)/10);
  float newMersf = ((clen-tailAmt)/10.0);
  float diff = newMersf - newMers;
  if (diff > 0) newMers++;

	// KJD Something is wrong with this.  It's allocating too much space. 
	//  cout<<"tailAmt="<<tailAmt<<"\tnewMersf="<<newMersf<<"\tdiff="<<diff<<"\tnewMers="<<newMers<<"\tnewSize=";
  //cout<<oldMerSize+newMers<<endl;

  mData.resize(oldMerSize+newMers);
  
  // Copy the bases to the Sequence3Bit vector. 
  for(idx = 0;idx<clen;idx++){
    if (mMerIdx < 10){
      (*this).Set(mBaseIdx,mMerIdx,cstring[idx]);
      mMerIdx++;
      mBaseSize++;
    }else{
      mMerIdx = 0;
      mBaseIdx++;
      (*this).Set(mBaseIdx,mMerIdx,cstring[idx]);
      mMerIdx++;
      mBaseSize++;
    }
  }
}

//----------------------------------------
// * Concatenate
//----------------------------------------
// Adds a Sequence3Bit classs sequence to this Sequence3Bit.   
//
void 
Sequence3Bit::Concatenate(Sequence3Bit &Seq)
{
  register int clen = Seq.size();
  register int oldMerSize = mData.size();
  register long idx;
  
  // Make space for it.
  int tailAmt = 10 - mMerIdx - 1;  
  int newMers = ((clen-tailAmt)/10);
  float newMersf = ((clen-tailAmt)/10.0);
  float diff = newMersf - newMers;
  if (diff > 0) newMers++;
  mData.resize(oldMerSize+newMers);
  
  // Copy the bases to the Sequence3Bit vector. 
  for(idx = 0;idx<clen;idx++){
    if (mMerIdx < 10){
      (*this).Set(mBaseIdx,mMerIdx,Seq[idx]);
      mMerIdx++;
      mBaseSize++;
    }else{
      mMerIdx = 0;
      mBaseIdx++;
      (*this).Set(mBaseIdx,mMerIdx,Seq[idx]);
      mMerIdx++;
      mBaseSize++;
    }
  }
}

//----------------------------------------
// * ConcatenateSequence
//----------------------------------------
// Adds a standard Sequence classs sequence to this Sequence3Bit.   
//
void 
Sequence3Bit::ConcatenateSequence(Sequence &Seq)
{
  register int clen = Seq.size();
  register int oldMerSize = mData.size();
  register long idx;
  
  // Make space for it.
  int tailAmt = 10 - mMerIdx - 1;  
  int newMers = ((clen-tailAmt)/10);
  float newMersf = ((clen-tailAmt)/10.0);
  float diff = newMersf - newMers;
  if (diff > 0) newMers++;
  mData.resize(oldMerSize+newMers);
  
  // Copy the bases to the Sequence3Bit vector. 
  for(idx = 0;idx<clen;idx++){
    if (mMerIdx < 10){
      (*this).Set(mBaseIdx,mMerIdx,Seq[idx]);
      mMerIdx++;
      mBaseSize++;
    }else{
      mMerIdx = 0;
      mBaseIdx++;
      (*this).Set(mBaseIdx,mMerIdx,Seq[idx]);
      mMerIdx++;
      mBaseSize++;
    }
  }
}

//----------------------------------------
// * ConcatenateCStringTrimLowQual
//----------------------------------------
// Adds a C-language string of ASCII bases to the Sequence3Bit, excluding the lower case
// (low quality soft trim) letters at the end. Assume that low-quality trim at begining of 
// read has been skipped before calling this function.  So, only looking for the beginning 
// of a new low-quality run.   Also assumes that lower case low quality trim is only in 
// a single block at beginning and end of sequence.  Any interiror lowercase bases will 
// cut off the sequence prematurely.  
//
void 
Sequence3Bit::ConcatenateCStringTrimLowQual(const char *cstring){
  register int clen = strlen(cstring);
  register int oldMerSize = mData.size();
  register long idx;
  
  // Make space for it.
  int tailAmt = 10 - mMerIdx - 1;  
  int newMers = ((clen-tailAmt)/10);
  float newMersf = ((clen-tailAmt)/10.0);
  float diff = newMersf - newMers;
 
  if (diff > 0) newMers++;
  // KJD BUG for some reason the calculations above fall short.  No harm to add one
  // for now. 
  newMers++;
  mData.resize(oldMerSize+newMers);
  
  // Copy the bases to the Sequence3Bit vector. 
  for(idx = 0;idx<clen;idx++){

    // Found a low quality base.  So, we're finished before end of this line. 
    if (cstring[idx] & 0x20){
      // If want to reclaim over estimated size needed. 
      //resize(oldMerSize+mBaseIdx)
      return;
    }
	 
    if (mMerIdx < 10){
      (*this).Set(mBaseIdx,mMerIdx,cstring[idx]);
      mMerIdx++;
      mBaseSize++;
    }else{
      mMerIdx = 0;
      mBaseIdx++;
      (*this).Set(mBaseIdx,mMerIdx,cstring[idx]);
      mMerIdx++;
      mBaseSize++;
    }
  }
}

//----------------------------------------
// * GetSubSequence
//----------------------------------------
void 
Sequence3Bit::GetSubSequence(Sequence3Bit &Seq,int start, int end){
  
  // Clear any previous sequence in the target.
  Seq.mData.clear();


  if (end < 0) return;

  // Truncate ends
  if (start < 0) start = 0;
  if (end >= this->size()) end = this->size()-1;
  
  long ssize = end - start + 1; 
  if (ssize > 0){
    //    Seq.resize(ssize);  // XXX Need a base level resize. 
    
    long subidx = 0;
    for(long idx = start;idx <= end;idx++){
      //  Seq[subidx++] = (*this)[idx];
      Seq += (*this)[idx];
    }

    string description = (*this).mDescription;
    // KJD Need to root these character arrays out. Suggest writing some
    // string utilities that handle int->string, float->string, etc. 
    char buf[500];
    sprintf(buf,"[%ld,%ld]",start,end);
    description = description+buf;
    Seq.mDescription = description;
    Seq.mName = (*this).mName;
  }
}

//----------------------------------------
// * GetNextReadTrimLowQual
//----------------------------------------
// Gets the next read from the given stream, omitting lower case (low quality) soft trim at 
// beginning and end of read.  This function works in tight conjunction with ConcatenateCStringTrimLowQual. 
// buffer[0] should be set to '1' for the first call of this function.  On subsequent calls, 
// buffer will contain the comment line for the next read in the file.  
//
bool 
Sequence3Bit::GetNextReadTrimLowQual(istream &fin,char *buffer,int bufSize)
{
  // Erase any existing sequence, just in case.  
  mData.clear();
    
  // Either we are at the first of file, indicated by '1', or we should
  // have a '>' for the next read.  Otherwise we are at the end of the file
  // or the file is corrupt or something. 
  if (buffer[0] == '1'){
    fin.getline(buffer,bufSize);
  }else{
    if (buffer[0] != '>') return(false); // Not at start of new read. 
  }
  
  // Extract name and description from line.  
  string nextDesc = buffer;
  int offset = 1; // To skip >
  // Everything up to first white space is name, after that is description. 
  mName = GetNextWSDelimitedString(nextDesc,offset,offset);  
  // It's possible there is no description, so check:
  if ((offset >= 0) && (offset < nextDesc.size())){
    mDescription = nextDesc.substr(offset,nextDesc.size()-offset);
  }else{
    mDescription = "";
  }

  // Find the first upper case (high quality) base in the sequence. 
  // KJD only works with standard ASCII and if sequence is guaranteed to be one of 
  // {a,c,g,t,n,x,A,C,G,T,N,X}  0x20 bit is set with lower case.  0x08 is set with n or x. 
  char c;
  do{
    fin.get(c);
  }while(!fin.eof() && (c & 0x20) && (c!= '>')); 

  // The entire read is low quality.  Now at next read.   
  if (c == '>'){	
    buffer[0] = '>';
    fin.getline(buffer+1,bufSize); // Get rest of comment line. 
    mBaseSize = 0;
    return(true);
  }

  bool bFirstTime = true;

  // Get the sequence until another > is encountered. 
  // KJD UGH!  Or until another low quality base is encountered!!!
  do{
    int bytesRead;
    if (bFirstTime){
      buffer[0] = c;
      fin.getline(buffer+1,bufSize);
      bytesRead = fin.gcount()+1;  // +1 is for c tacked onto the beginning. 
      bFirstTime = false;
    }else{
      fin.getline(buffer,bufSize);
      bytesRead = fin.gcount();
    }
    int bsize = strlen(buffer);

    // This means that we have a null before the newline, so I abort this
    // probably corrupt read.   I found a few of these and it was a pain to 
    // have my program choke 1 million reads into some operation because one read was bad. 
    //if ((bsize < (bytesRead-1)) && (bsize !=0) && (bytesRead != 0)) return(false);
    // KJD I'm going to assume that this doesn't happen in our current data.  
	 
    if (buffer[0] != '>'){
      // Sadly, many FASTA files have empty lines or whitespace at the 
      // end of lines.  This is effectively filtered here by adding a null
      // terminator in place of all whitespace. 
		
      register int sEnd;
      if (bsize > 0) sEnd = bsize-1;
      else sEnd = 0;
      while((buffer[sEnd] == '\t') || 
	    (buffer[sEnd] == ' ') ||
	    (buffer[sEnd]=='\n') && (sEnd > 0)){
	cout<<"bsize="<<bsize<<endl;
	cout<<"Buffer="<<buffer<<endl;
	cout<<"Consuming white space: "<<hex<<(int)buffer[sEnd]<<" sEnd="<<sEnd<<endl;
	buffer[sEnd--] = '\0';
      }
      // Copy the characters of the line to the sequence vector. 		
      ConcatenateCStringTrimLowQual(buffer);			            
    }
  }while((!fin.eof()) && (buffer[0] != '>'));

  return(true);
}

//----------------------------------------
// * GetNextRead
//----------------------------------------
// Gets the next read from a file. 
// buffer[0] should be set to '1' for the first call of this function.  On subsequent calls, 
// buffer will contain the comment line for the next read in the file.  
// 
bool 
Sequence3Bit::GetNextRead(istream &fin,char *buffer,int bufSize)
{
  
  // Erase any existing sequence, just in case.  
  mData.clear();
  
  // Must be first line  
  if (buffer[0] == '1'){
    fin.getline(buffer,bufSize);
  }
  
  // Extract name and description from line.  
  string nextDesc = buffer;
  int offset = 1; // To skip >
  // Everything up to first white space is name, after that is description. 
  mName = GetNextWSDelimitedString(nextDesc,offset,offset);  
  // It's possible there is no description, so check:
  if ((offset >= 0) && (offset < nextDesc.size())){
    mDescription = nextDesc.substr(offset,nextDesc.size()-offset);
  }else{
    mDescription = "";
  }
  
  // Get the sequence until another > is encountered. 
  do{
    fin.getline(buffer,bufSize);
    string buf = buffer;
    int bytesRead = fin.gcount();
    int bsize = strlen(buffer);
    // This means that we have a null before the newline, so I abort this
    // probably corrupt read. 
    if ((bsize < (bytesRead-1)) && (bsize !=0) && (bytesRead != 0)) return(false);
	 
    if (buffer[0] != '>'){
      // Sadly, many FASTA files have empty lines or whitespace at the 
      // end of lines.  This is effectively filtered here by adding a null
      // terminator in place of all whitespace. 
		
      register int sEnd;
      if (bsize > 0) sEnd = bsize-1;
      else sEnd = 0;
      while((buffer[sEnd] == '\t') || 
	    (buffer[sEnd] == ' ') ||
	    (buffer[sEnd]=='\n') && (sEnd > 0)){
	cout<<"\nConsuming white space: "<<hex<<(int)buffer[sEnd]<<" sEnd="<<sEnd<<dec<<endl;
	cout<<"Buffer="<<buffer<<endl;
	buffer[sEnd--] = '\0';
      }

      // Copy the characters of the line to the sequence vector. 
		
      ConcatenateCString(buffer);			            
    }
  }while((!fin.eof()) && (buffer[0] != '>'));
  return(true);
}



